      program btu751
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      parameter (np1=4,mp1=np1+1,np2=7,mp2=np2+1)
      parameter (ftol=1.d-10,nrep=9)
c
      external famoeb,gamoeb
      character*6 name(np2+1)
      character*24 fdate
      dimension x1(np1),p1(mp1,np1),y1(mp1),x0(np1),x10(np2-np1)
      dimension x2(np2),p2(mp2,np2),y2(mp2)
      dimension is(nplay,ndec)
      dimension est1(np1),est2(np2),xx1(np1),xx2(np2)
      dimension bllf2(nrep),boot2(np2,nrep),bratio(nrep),
     &          btmp(nrep),bs2(np2),bss2(np2)
      common /mats/is
      common /est/est1
c
      data name/'BETA1','PZ1','SB1','SZ1','SB2','SZ2','A1','A2'/
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
      open(unit=10,file='btu751.dat')
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c-----------------------------------------------------------------
      write(*,'('' ***** bt9u51 ******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/" All 2011 Sessions"//)')
c
         call fgrid
         call fch
c
c Set parameters
         x0(1) = dlog(1.d0/0.99d0 - 1.d0)
         x0(2) = dlog(0.25d0/0.13d0 - 1.d0) 
         x0(3) = dlog(1.d0/0.21d0 - 1.d0)
         x0(4) = dlog(1.d0/0.14d0 - 1.d0)
c
         delta=0.1d0
         iter=3000  
         call init(x0,p1,y1,famoeb,delta)
         call amoeba(p1,y1,mp1,np1,np1,ftol,famoeb,iter)
         write(*,'("iterations = ",i5)')iter
         call avg(p1,x1,np1,mp1)
         z1 = -famoeb(x1)
         call m2h(x1,est1,np1)
         write(*,'("est = ",4g13.5)')(est1(k),k=1,np1)
         write(*,'(/"Initial LL = ",g14.6)')z1
c         write(10,910)0,0,z1,(est1(k),k=1,np1)
c         stop
c 
         call savephi(est1)
c
         x0(1) = dlog(1.d0/0.99d0 - 1.d0)
         x0(2) = dlog(0.25d0/0.11d0 - 1.d0) 
         x0(3) = dlog(1.d0/0.26d0 - 1.d0)
         x0(4) = dlog(1.d0/0.11d0 - 1.d0)
         x10(1) = dlog(1.d0/0.01d0 - 1.d0)
         x10(2) = dlog(1.d1/9.0d0 - 1.d0)
         x10(3) = dlog(1.d0/0.29d0 - 1.d0)
c
         delta=0.25d0
         iter=3000
         call jnit(x0,x10,p2,y2,gamoeb,delta)
         call amoeba(p2,y2,mp2,np2,np2,ftol,gamoeb,iter)
         call avg(p2,x2,np2,mp2)
         z2 = -gamoeb(x2)
         br0 = z2 - z1
         call m2h2(x2,xx2,np2)
         e2 = 1.d0-xx2(np2)
         write(*,'("iterations = ",i5)')iter
         write(*,'("est2 = ",8g13.5)')(xx2(k),k=1,np2),e2
         write(*,'(/"Initial LL = ",2g14.6)')z2,br0
c         stop
c
c begin interations
c
c
      dseed = 90001.d0
      do 500 it=1,nrep
   25    xseed = dseed
         call data(it,dseed)
         call fch
c
c uni model
c
         x0(1) = dlog(1.d0/0.99d0 - 1.d0)
         x0(2) = dlog(0.25d0/0.13d0 - 1.d0) 
         x0(3) = dlog(1.d0/0.21d0 - 1.d0)
         x0(4) = dlog(1.d0/0.14d0 - 1.d0)
         delta=0.5d0
         iter=3000
         call init(x0,p1,y1,famoeb,delta)
         call amoeba(p1,y1,mp1,np1,np1,ftol,famoeb,iter)
         if (iter.ge.3000)  then
           write(*,'(5x,"irep =",2x,i5)')it
           go to 25
         endif
         call avg(p1,x1,np1,mp1)
         z1 = -famoeb(x1)
c         call m2h(x1,xx1,np1)
c
c ramix9 model
         x0(1) = dlog(1.d0/0.99d0 - 1.d0)
         x0(2) = dlog(0.25d0/0.11d0 - 1.d0) 
         x0(3) = dlog(1.d0/0.26d0 - 1.d0)
         x0(4) = dlog(1.d0/0.11d0 - 1.d0) 

         x10(1) = dlog(1.d0/0.01d0 - 1.d0)
         x10(2) = dlog(1.d1/9.9d0 - 1.d0)
         x10(3) = dlog(1.d0/0.29d0 - 1.d0)
         iflag=0
   22    continue
         delta=0.5d0
         iter=5000
         call jnit(x0,x10,p2,y2,gamoeb,delta)
         call amoeba(p2,y2,mp2,np2,np2,ftol,gamoeb,iter)
         if (iter.ge.5000)  then
           write(*,'(5x,"irep =",2x,i5)')it
c           do 21 i=1,mp2
c             write(*,'(10g12.4)')(p2(i,j),j=1,np2),y2(i)
c   21      continue
           go to 25
         endif
         call avg(p2,x2,np2,mp2)
         bllf2(it) = -gamoeb(x2)
         bratio(it) = bllf2(it) - z1
c         if ((bratio(it) .lt. -1.d0).and.(iflag.eq.0)) then
c           write(*,'(5x,i4," bratio < -1")')it
c           x0(1) = dlog(1.d0/0.99d0 - 1.d0)
c           x0(2) = dlog(0.25d0/0.01d0 - 1.d0) 
c           x0(3) = dlog(1.d0/0.16d0 - 1.d0)
c           x0(4) = dlog(1.d0/0.27d0 - 1.d0)
c           x10(3) = dlog(1.d0/0.01d0 - 1.d0)
c           iflag=1
c           goto 22
c         endif
         call m2h2(x2,xx2,np2)
         do 30 k=1,np2
            boot2(k,it) = xx2(k)
   30    continue
         write(10,910)it,xseed,bllf2(it),(boot2(k,it),k=1,np2)
         write(*,911)it,z1,bllf2(it),bratio(it)
  500 continue
  910   format(i4,2x,f20.1,5(2x,g14.6))
  911   format(i3,2x,3(2x,g14.6))
c
c --- process the bootstrap results ---
c
      call sort(nrep,bratio)
c      do 510 it=1,nrep
c         write(*,'(i4,x,2g13.5)')it,bratio(it)
c  510 continue
      nhi = nrep - 2.d0*t1
      write(*,'(//"The 10% critical llratio is ",g12.4)')bratio(nhi)
      do 520 i=1,nrep
         i0 = nrep - i + 1
         if (bratio(i0).lt.br0) goto 521
  520 continue
  521 pv = float(i-1)/float(nrep)
      write(*,'(/"The p-value of ",f7.3," is ",g12.4)')br0,pv
c
      end
c***********************************************************************
      subroutine m2h(x,xx,np)
      implicit real*8 (a-h,o-z)
      parameter (ngrid=51)
      dimension x(np),xx(np)
c
      s0=0.05d0
      xx(1) = 1.d0/(1.d0+dexp(x(1)))
      xx(2) = 0.25d0/(1.d0+dexp(x(2)))+0.5d0
      xx(3) = 1.d0/(1.d0+dexp(x(3)))+s0
      xx(4) = 1.d0/(1.d0+dexp(x(4)))+s0
c
      return
      end
c***********************************************************************
      subroutine m2h2(x,xx,np)
      implicit real*8 (a-h,o-z)
      parameter (ngrid=51)
      dimension x(np),xx(np)
c
      s0=0.05d0
      xx(1) = 1.d0/(1.d0+dexp(x(1)))
      xx(2) = 0.25d0/(1.d0+dexp(x(2)))+0.5d0
      xx(3) = 1.d0/(1.d0+dexp(x(3)))+s0
      xx(4) = 1.d0/(1.d0+dexp(x(4)))+s0
      xx(5) = 1.d0/(1.d0+dexp(x(5)))+s0
      xx(6) = 9.95d0/(1.d0+dexp(x(6)))+s0
      xx(7) = 1.d0/(1.d0+dexp(x(7)))
c
      return
      end
c***********************************************************************
      double precision function famoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=4,ngrid=51)
      dimension x(np)
c
      s0=0.05d0
      beta = 1.d0/(1.d0+dexp(x(1)))
      pz = 0.25d0/(1.d0+dexp(x(2)))+0.5d0
      sb=1.d0/(1.d0+dexp(x(3)))+s0
      sz=1.d0/(1.d0+dexp(x(4)))+s0
c
      call fam(beta,pz,sb,sz,pt)

      famoeb = -pt
c
      return
      end
c***********************************************************************
      subroutine fam(beta,pz,sb,sz,pt)
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pch(nplay,ngrid,ngrid),phi(ngrid,ngrid)
      common /pchm/pch

      call fprob(beta,pz,sb,sz,phi)
c
      pt = 0.d0
      do 50 id=1,nplay
         pid = 0.d0
         do 20 j=1,ngrid
            do 15 k=1,ngrid
               pid = pid + phi(j,k)*pch(id,j,k)
   15       continue
   20    continue
         pid = dlog(pid)
         pt = pt + pid
   50 continue
      return
      end
c***********************************************************************
      subroutine fgrid
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension w0(ngrid),w(ngrid,ngrid)
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      common /grid/pxm,w
c
      do 5 i=1,ngrid
         w0(i) = 1.d0
         if((i.eq.1).or.(i.eq.ngrid)) w0(i) = 0.5d0
    5 continue
c
      dz=1.d0/float(ngrid-1)
      do 30 j=1,ngrid
         beta = dz*(j-1)
         do 20 k=1,ngrid
            if (k.lt.ngrid) then
               pz = 0.5d0 + 0.5d0*dz*(k-1)
            else
               pz = 0.999d0
            endif
            gam = 2.d0*dlog(pz/(1.d0-pz))
            call fpx(gam,beta,px)
            do 10 m=1,ndec
               pxm(j,k,m,1) = px(m,1)
               pxm(j,k,m,2) = px(m,2)
   10       continue
            w(j,k) = w0(j)*w0(k)
   20    continue
   30 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fch
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension pch(nplay,ngrid,ngrid)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /mats/is
      common /pchm/pch

      do 50 id=1,nplay
         do 30 j=1,ngrid
            do 25 k=1,ngrid
               pin = 1.d0
               do 20 m=1,ndec
                  ix = is(id,m)
                  pin = pin*pxm(j,k,m,ix)
   20          continue
               pch(id,j,k) = pin
   25       continue
   30    continue
   50 continue
      return
      end
c***********************************************************************
      subroutine fprob(beta,pz,sb,sz,phi)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension phi(ngrid,ngrid)
      common /grid/pxm,w
c
      pt = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 20 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - beta)/sb
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz)/sz
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi(j,k) = w(j,k)*pg*pb
            pt = pt + phi(j,k)
   10    continue
   20 continue
      do 30 j=1,ngrid
         do 25 k=1,ngrid
            phi(j,k) = phi(j,k)/pt
   25    continue
c      write(*,'(11f7.3)')(phi(j,k),k=1,ngrid)
   30 continue
c      stop
      return
      end
c***********************************************************************
      subroutine savephi(x)
      implicit real*8 (a-h,o-z)
      parameter (np=4,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension x(np),gl(npts)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      common /grid/pxm,w
      common /glm/gl
c
      beta = x(1)
      pz = x(2)
      sb = x(3)
      sz = x(4)
c      write(*,'(4g12.4)')beta,pz,sb,sz
c
      gt = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 20 j=1,ngrid
         j0 = ngrid*(j-1)
         bj = dz*(j-1)
         pb = (bj - beta)/sb
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz)/sz
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            gl(j0+k) = w(j,k)*pg*pb
            gt = gt + gl(j0+k)
   10    continue
   20 continue
c
      gl(1) = gl(1)/gt
      do 30 i=2,npts
         gl(i) = gl(i-1) + gl(i)/gt
   30 continue
c      write(*,'(21f8.4)')(gl(i),i=1,npts)
c      stop
      return
      end
c***********************************************************************
      double precision function gamoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=9,ngrid=51)
      dimension x(np)
c
      s0=0.05d0
      b1  = 1.d0/(1.d0+dexp(x(1)))
      pz1 = 0.25d0/(1.d0+dexp(x(2)))+0.5d0
      sb1 = 1.d0/(1.d0+dexp(x(3)))+s0
      sz1 = 1.d0/(1.d0+dexp(x(4)))+s0
      b2  = 1.d0
      pz2 = 0.999d0
      sb2 = 1.d0/(1.d0+dexp(x(5)))+s0
      sz2 = 9.95d0/(1.d0+dexp(x(6)))+s0
      a1  = 1.d0/(1.d0+dexp(x(7)))
c
      call gam(b1,pz1,sb1,sz1,b2,pz2,sb2,sz2,a1,pt)

      gamoeb = -pt
c
      return
      end
c***********************************************************************
      subroutine gam(b1,pz1,sb1,sz1,b2,pz2,sb2,sz2,a2,pt)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pch(nplay,ngrid,ngrid),phi1(ngrid,ngrid),
     &          phi2(ngrid,ngrid)
      common /pchm/pch

      call gprob(b1,pz1,sb1,sz1,b2,pz2,sb2,sz2,phi1,phi2)
      a1 = 1.d0-a2
c
      pt = 0.d0
      do 50 id=1,nplay
         pid = 0.d0
         do 20 j=1,ngrid
            do 15 k=1,ngrid
               phi = a1*phi1(j,k) + a2*phi2(j,k)
               pid = pid + phi*pch(id,j,k)
   15       continue
   20    continue
         if (pid .gt. 1.d-307) then
            pid = dlog(pid)
         else
            pid = -706.d0
         endif
         pt = pt + pid
   50 continue
      return
      end
c***********************************************************************
      subroutine gprob(b1,pz1,sb1,sz1,b2,pz2,sb2,sz2,phi1,phi2)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension phi1(ngrid,ngrid),phi2(ngrid,ngrid)
      common /grid/pxm,w
c
      dz = 1.d0/float(ngrid-1)
      pt1 = 0.d0
      do 20 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - b1)/sb1
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz1)/sz1
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi1(j,k) = w(j,k)*pg*pb
            pt1 = pt1 + phi1(j,k)
   10    continue
   20 continue
      do 30 j=1,ngrid
         do 25 k=1,ngrid
            phi1(j,k) = phi1(j,k)/pt1
   25    continue
c      write(*,'(11f7.3)')(phi1(j,k),k=1,ngrid)
   30 continue
c
      pt2 = 0.d0
      do 50 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - b2)/sb2
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 40 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz2)/sz2
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi2(j,k) = w(j,k)*pg*pb
            pt2 = pt2 + phi2(j,k)
   40    continue
   50 continue
      do 60 j=1,ngrid
         do 55 k=1,ngrid
            phi2(j,k) = phi2(j,k)/pt2
   55    continue
c      write(*,'(11f7.3)')(phi2(j,k),k=1,ngrid)
   60 continue
c
      return
      end
c***************************************************************
      subroutine data(it,dseed)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension gl(npts)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /glm/gl
      common /mats/is
c
      do 50 id=1,nplay
         r = ggubfs(dseed)
         do 15 i=1,npts
            if (gl(i) .gt. r) goto 16
   15    continue
c
   16    j = (i-1)/ngrid + 1
         k = i - ngrid*(j-1)
c
         do 30 m=1,ndec
            x = ggubfs(dseed)
            p1 = pxm(j,k,m,1)
            if (x .le. p1) then 
               is(id,m)=1
            else
               is(id,m)=2
            endif
   30    continue
c      write(*,'(13i3)')id,(is(id,k),k=1,ndec)
   50 continue
      return
      end
c
c***********************************************************************
      subroutine init(xint,p,y,famoeb,delta)
      implicit real*8 (a-h,o-z) 
      parameter (np=4,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-famoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
c      delta=0.1d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c***********************************************************************
      subroutine jnit(x1,x2,p,y,gamoeb,delta)
      implicit real*8 (a-h,o-z) 
      parameter (kp=4,np=7,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension x1(kp),x2(np-kp),xx(np),xxx(np)
c 

      do 1 k=1,kp-1
         xxx(k) = x1(k)
         xxx(k+4) = x2(k)
    1 continue
      xxx(kp) = x1(kp)
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-gamoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=gamoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end
c**************************************************************
      subroutine sort(n,ra)
      implicit real*8 (a-h,o-z)
      dimension ra(n)
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
        l=l-1
        rra=ra(l)
        else
        rra=ra(ir)
        ra(ir)=ra(1)
        ir=ir-1
        if(ir.eq.1)then
        ra(1)=rra
        return
        endif
        endif
        i=l  
        j=l+l
20      if(j.le.ir)then
        if(j.lt.ir)then
        if(ra(j).lt.ra(j+1))j=j+1
        endif
        if(rra.lt.ra(j))then
        ra(i)=ra(j)
        i=j
        j=j+j
        else
        j=ir+1
        endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end      
c========================================================================
      double precision function ggubfs (dseed)
      double precision   dseed
      double precision   d2p31m,d2p31
      data              d2p31m/2147483647.d0/
      data              d2p31 /2147483648.d0/
      dseed = dmod(16807.d0*dseed,d2p31m)
      ggubfs = dseed / d2p31
      return
      end

